Introduction

The experimental design is described in the preregistration here: https://osf.io/m9xep (@Daniela, you can skip the section Lab Equipment)

In this markdown, we will be modelling the causal relation between effort and correction to confirm/reject our hypothesis.

These are:

H1: In corrections, people will enhance some effort-related kinematic and/or acoustic features of their behaviour

H2: The enhancement will depend on similarity of the guesser’s answer and the original meaning. More similar answer will require/result in smaller enhancement (but still enhancement) than less similar answer.

Setting up the environment

library(here)
## here() starts at E:/FLESH_ContinuousBodilyEffort
library(dagitty) # for dag
library(dplyr) # for data-wrangling
## 
## Attache Paket: 'dplyr'
## Die folgenden Objekte sind maskiert von 'package:stats':
## 
##     filter, lag
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lme4)  # for linear mixed-effects models
## Lade nötiges Paket: Matrix
library(tidyr)  # for reshaping data (if needed)
## 
## Attache Paket: 'tidyr'
## Die folgenden Objekte sind maskiert von 'package:Matrix':
## 
##     expand, pack, unpack
library(ggplot2)

DAG

To make sure we know what variables should be in the model and why, we draw a DAG that illustrates our assumptions

Our relationship of interest is the causal effect of communicative attempt on effort.

Other assumptions include:

Personality traits (measured with Big5) will influence effort (e.g., people are more willing to put more effort if they are open-minded) and also communicative attempt (e.g., more extroverted people are better in this game, therefore they need less attempts)

Familiarity with guessing partner influences effort (ditto) as well as communicative attempt (e.g., friends are better in this game than strangers)

Similarly, participant will also directly infleunce effort and commAtt, because personality traits might not be capturing all variation

Expressibility of concepts is going to influence effort (e.g., more expressible concepts allow more enhancement - but could be also other direction) as well as CommAtt (e.g., more expressible concepts are more readily guessed and don’t need more attempts)

Similarly, concept will also directly influence effort and commAtt, because expressibility might not be capturing all variation

Modality (uni vs. multi) will influence the value of effort. We assume that in unimodal condition, the feature does not need to account for synergic relations with the other modality, and may carry the whole signal. In multimodal condition, however, there may be synergic relations that moderate the full expressiveness of this feature

Lastly, trial number (characterising how for one is in the experiment) could be hinting on learning processes through out the experiment, or - the other direction - on increasing fatigue

dag <- dagitty('dag {
Big5 [adjusted,pos="-0.823,0.657"]
CommAtt [exposure,pos="-1.033,0.028"]
Conc [adjusted,pos="-1.136,-0.848"]
Eff [outcome,pos="-0.102,0.025"]
Expr [adjusted,pos="-0.758,-0.850"]
Fam [adjusted,pos="-0.379,0.663"]
Mod_cat [adjusted,pos="-0.374,-0.850"]
Pcn [adjusted,pos="-0.589,1.214"]
TrNum [adjusted,pos="-1.686,-0.859"]
Big5 -> CommAtt     
Big5 -> Eff
CommAtt -> Eff
Conc -> Expr
Expr -> CommAtt
Expr -> Eff
Fam -> CommAtt
Fam -> Eff
Mod_bin -> Eff
Mod_cat -> Expr
Pcn -> Big5
Pcn -> CommAtt 
Pcn -> Eff
Pcn -> Fam
TrNum -> CommAtt
TrNum -> Eff
Conc -> CommAtt
Conc -> Eff
}')

plot(dag)
## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.

Synthetic data

We will now create synth data that will copy the relations we assume in our DAG. Assigning concrete coefficients will also help us to test our model - we should find exactly those causalities that we had in mind when creating these data

# Set seed for reproducibility
set.seed(0209)

# Define participants, total unique concepts, and modalities
n_participants <- 120
n_total_concepts <- 21  # Total unique concepts
n_concepts_per_participant <- 21  # Each participant works with 21 concepts
n_modalities <- 3  # gesture, vocal, combined

# Generate participant IDs
participants <- 1:n_participants

# Simulate Big5 personality traits (standardized between 0 and 1) and Familiarity (between 0 and 1) for participants
Big5 <- runif(n_participants, min = 0, max = 1)  # Continuous values between 0 and 1
Familiarity <- runif(n_participants, min = 0, max = 1)  # Continuous values between 0 and 1

# Create a matrix to hold expressibility values for each concept in each modality
expressibility_matrix <- matrix(runif(n_total_concepts * n_modalities, min = 0, max = 1), nrow = n_total_concepts, ncol = n_modalities)

# Randomly sample 21 unique concepts for each participant
final_data <- data.frame()

# Define a function to assign CommAtt and Eff for a single participant
simulate_participant <- function(participant_id) {
  # Randomly sample 21 unique concepts from the total pool of 84
  selected_concepts <- sample(1:n_total_concepts, n_concepts_per_participant)
  
  participant_data <- data.frame()
  trial_number <- 1  # Initialize trial number
  
  for (concept_id in selected_concepts) {
    # Randomly determine the modality for the concept
    modality <- sample(c("gesture", "vocal", "combined"), 1)
    
    # Calculate expressibility based on modality
    expressibility_score <- ifelse(modality == "vocal", expressibility_matrix[concept_id, 1] * 0.6, 
                                    ifelse(modality == "gesture", expressibility_matrix[concept_id, 2], 
                                           expressibility_matrix[concept_id, 3] * 1.2))
    
    # Determine Communicative Attempts based solely on expressibility, familiarity, and Big5
    base_prob <- c(0.33, 0.33, 0.33)  # Equal chance for 1, 2, or 3 attempts
    
    # Modify probabilities based on familiarity, Big5, and expressibility
    adjusted_prob <- base_prob * c(1 - Familiarity[participant_id], # 3 times for each
                                    1 - Familiarity[participant_id],
                                    1 - Familiarity[participant_id]) * 
                     c(1 - Big5[participant_id],
                       1 - Big5[participant_id],
                       1 - Big5[participant_id]) * 
                     c(1 - expressibility_score,
                       1 - expressibility_score,
                       1 - expressibility_score)
    
    # Normalize the adjusted probabilities
    adjusted_prob <- adjusted_prob / sum(adjusted_prob)
    
    # Sample the number of communicative attempts based on adjusted probabilities
    n_attempts <- sample(1:3, 1, prob = adjusted_prob)
    
    # Loop through the number of attempts and increment CommAtt correctly
    for (attempt in 1:n_attempts) {
      # Calculate Eff for the first attempt
      if (attempt == 1) {
        Eff <- 1.15 * Big5[participant_id] + 
               1.10 * Familiarity[participant_id] + 
               1.20 * expressibility_score + 
               rnorm(1, mean = 1, sd = 0.5)
        
        # Adjust Eff based on modality
        if (modality == "combined") {
          
          Eff <- Eff * 0.7  # Slight moderation for combined modality
        }
      }
      
      # Adjust Eff for subsequent attempts
      if (attempt == 2) {
        Eff <- 1.15 * Big5[participant_id] + 
               1.10 * Familiarity[participant_id] + 
               1.20 * expressibility_score + 
               rnorm(1, mean = 1, sd = 0.5)
        Eff <- Eff * 1.50  # Multiply effort by 1.50 for the second attempt
      } else if (attempt == 3) {
        Eff <- 1.15 * Big5[participant_id] + 
               1.10 * Familiarity[participant_id] + 
               1.20 * expressibility_score + 
               rnorm(1, mean = 1, sd = 0.5)
        Eff <- Eff * 0.70  # Multiply effort by 0.70 for the third attempt
      }
      
      # Create row for each attempt
      participant_data <- rbind(participant_data, data.frame(
        Participant = participant_id,
        Concept = concept_id,
        Modality = modality,
        Big5 = Big5[participant_id],
        Familiarity = Familiarity[participant_id],
        Expressibility = expressibility_score,
        CommAtt = attempt,  # Correctly set the attempt number
        Eff = Eff,
        TrialNumber = trial_number  # Set trial number for this attempt
      ))
      
      # Increment the trial number after each attempt
      trial_number <- trial_number + 1
    }
  }
  
  return(participant_data)
}

# Simulate data for all participants
for (i in participants) {
  final_data <- rbind(final_data, simulate_participant(i))
}

# Preview the first few rows of the final data
head(final_data)
##   Participant Concept Modality      Big5 Familiarity Expressibility CommAtt
## 1           1       3 combined 0.9599493   0.9219529     0.07567109       1
## 2           1       3 combined 0.9599493   0.9219529     0.07567109       2
## 3           1       3 combined 0.9599493   0.9219529     0.07567109       3
## 4           1       9    vocal 0.9599493   0.9219529     0.37757391       1
## 5           1       9    vocal 0.9599493   0.9219529     0.37757391       2
## 6           1      21  gesture 0.9599493   0.9219529     0.56918606       1
##        Eff TrialNumber
## 1 2.624256           1
## 2 2.715052           2
## 3 2.650435           3
## 4 3.498290           4
## 5 5.902036           5
## 6 3.246798           6

So now we have synthetic data where we exactly defined (using coefficients) what is the relationship between certain variables

  1. CommAtt -> Eff

The effort for second attempt is multiplied by 1.25 (Beta = 1.25) The effort for third attempts by 0.9 (ie decreases, Beta = 0.9)

  1. Fam -> Eff & CommAtt

Beta = 1.10 for Eff Negative effect on CommAtt

  1. Big5 -> Eff & CommAtt

Beta = 1.15 for Eff Negative effect on CommAtt

  1. Expr -> Eff & CommAtt

Beta = 1.20 for Eff Negative effect on CommAtt

  1. TrNum -> Eff & CommAtt

Beta = could be both neg (fatigue) or positive (learning)

Exploring structure (DP)

nrow(final_data)
## [1] 5055
final_data |> 
  janitor::tabyl(Participant, Concept)
##  Participant 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
##            1 3 2 3 3 3 1 1 1 2  1  2  2  3  1  1  2  1  1  3  2  2
##            2 2 1 2 3 3 3 2 1 2  1  1  2  1  1  3  3  2  2  2  3  3
##            3 3 3 3 2 1 3 2 1 3  1  2  3  1  3  3  3  2  3  1  1  3
##            4 3 3 3 1 2 3 2 2 1  3  1  3  1  1  2  3  1  2  2  2  2
##            5 2 1 2 1 2 3 3 1 2  1  3  2  1  3  3  2  1  2  2  2  3
##            6 2 2 3 3 3 3 1 2 3  2  1  3  2  1  2  1  2  3  2  1  3
##            7 2 3 3 2 2 3 3 3 3  1  3  3  2  3  3  2  2  2  3  3  3
##            8 1 2 2 1 2 2 1 3 2  2  2  3  1  1  3  3  2  3  1  3  1
##            9 1 1 1 2 3 1 3 2 1  1  2  1  2  3  3  2  1  3  3  3  1
##           10 2 1 3 2 3 2 1 1 2  2  2  3  3  1  1  3  3  3  1  3  2
##           11 3 2 3 3 3 3 1 2 3  3  1  1  3  2  3  2  2  3  3  1  2
##           12 2 3 3 2 3 2 3 3 2  1  3  1  1  3  3  1  2  1  1  1  3
##           13 1 2 2 1 3 1 3 1 2  2  3  3  3  1  3  3  2  1  2  1  1
##           14 3 2 2 2 1 2 3 2 2  3  1  2  3  2  1  1  3  3  2  1  2
##           15 3 1 2 1 2 2 2 3 1  3  2  1  1  2  2  3  3  1  3  1  3
##           16 2 3 1 2 2 3 1 3 3  3  2  3  2  2  3  2  3  1  2  2  3
##           17 3 2 3 1 3 2 3 3 3  2  2  3  3  1  3  2  3  2  1  1  1
##           18 3 3 2 1 2 3 2 1 1  1  1  1  3  2  1  1  3  3  2  2  2
##           19 3 1 1 1 1 2 2 3 3  1  2  1  2  2  1  1  2  2  3  3  3
##           20 1 1 3 1 1 1 2 2 3  3  3  2  3  2  2  3  3  3  2  3  3
##           21 3 1 1 3 2 3 1 1 3  3  2  3  2  3  1  1  2  1  3  1  2
##           22 1 3 1 2 1 3 3 1 3  2  3  1  1  1  1  3  1  1  2  1  3
##           23 3 3 2 3 2 1 1 1 2  1  2  3  2  3  2  3  3  1  2  1  1
##           24 2 1 3 1 2 2 3 1 3  3  1  3  3  2  2  1  1  3  2  3  2
##           25 2 3 2 3 3 3 2 2 1  3  3  3  1  3  2  3  3  2  1  3  2
##           26 3 1 3 3 2 3 3 1 2  1  2  3  1  1  2  1  3  3  1  2  1
##           27 1 1 3 1 3 3 1 2 3  2  3  3  1  2  1  2  3  2  3  1  3
##           28 2 2 1 1 2 3 1 1 2  3  3  1  3  2  2  3  1  1  2  1  1
##           29 2 3 2 1 1 3 3 1 3  3  1  2  1  2  2  2  3  1  2  1  1
##           30 2 3 2 2 3 2 3 1 1  3  2  1  2  2  2  2  2  1  1  3  2
##           31 3 3 2 3 3 2 1 1 2  3  1  1  1  3  3  3  2  1  3  3  2
##           32 1 3 3 3 2 1 2 3 2  3  3  2  1  2  3  1  2  1  3  3  1
##           33 2 2 1 1 1 1 2 2 3  3  2  2  3  2  1  3  2  3  1  3  2
##           34 3 1 3 1 3 3 2 3 1  2  3  3  2  3  1  1  3  1  3  2  1
##           35 2 3 1 3 3 2 2 1 1  1  2  2  1  1  3  1  2  2  2  3  3
##           36 3 2 1 2 2 2 3 2 2  1  2  2  1  3  3  3  3  3  1  3  3
##           37 2 1 2 2 2 1 3 1 2  1  1  1  3  3  2  2  1  3  3  1  3
##           38 2 2 3 1 3 2 1 2 1  2  3  1  1  2  1  1  3  1  2  3  3
##           39 1 2 3 1 2 1 1 2 1  3  1  1  1  2  1  1  1  2  2  3  3
##           40 1 1 3 2 3 2 1 3 1  2  1  3  3  3  3  3  3  2  3  1  2
##           41 2 2 1 1 1 2 2 3 2  3  1  3  1  2  3  3  1  3  3  2  2
##           42 3 2 1 3 3 1 3 2 2  3  1  1  2  1  1  3  1  1  1  2  3
##           43 2 3 2 3 3 2 3 2 2  2  3  1  1  2  2  3  2  1  1  3  2
##           44 3 1 1 1 1 1 3 1 2  1  3  3  1  2  2  1  3  3  3  2  3
##           45 2 1 1 3 2 1 1 1 3  2  2  3  1  3  1  2  1  2  2  1  1
##           46 2 3 3 3 2 2 2 1 2  3  2  3  3  1  2  1  3  2  3  1  1
##           47 1 3 1 2 1 3 1 2 3  3  3  3  2  3  3  2  3  3  3  2  2
##           48 1 1 1 3 1 1 2 1 3  2  1  2  3  3  1  1  3  3  3  2  2
##           49 1 2 3 2 3 3 1 2 2  2  1  1  2  2  2  3  2  3  2  3  2
##           50 1 1 1 2 2 1 1 1 1  2  2  1  1  1  1  3  1  2  3  2  3
##           51 1 2 3 2 3 1 3 3 2  3  1  2  2  2  3  1  2  2  1  2  1
##           52 3 1 3 3 2 2 2 3 3  2  2  2  2  2  3  1  2  3  3  1  2
##           53 1 1 2 3 2 1 2 3 2  2  3  1  1  3  3  3  1  1  1  1  2
##           54 2 3 3 2 1 1 2 1 2  1  2  1  2  1  1  1  1  3  2  3  3
##           55 1 2 2 3 1 2 1 2 1  3  1  2  3  3  1  2  1  3  2  3  3
##           56 3 1 2 2 1 1 2 3 1  1  3  1  2  2  1  3  3  1  2  2  2
##           57 3 3 1 1 3 2 3 3 1  2  1  1  3  1  3  3  3  1  1  2  3
##           58 3 2 3 1 1 3 3 1 1  2  3  3  1  2  1  1  2  2  3  3  1
##           59 2 1 2 1 1 2 2 3 3  2  1  2  3  1  3  2  3  1  2  2  1
##           60 2 3 2 2 3 2 3 1 2  2  3  1  3  3  1  1  1  1  1  1  1
##           61 1 1 2 2 1 2 2 1 3  3  1  1  3  2  1  1  3  3  1  1  2
##           62 1 1 2 1 1 2 1 3 1  1  2  1  3  3  1  2  2  3  2  3  3
##           63 2 3 3 2 2 2 1 3 1  2  1  3  1  3  3  3  1  1  2  2  2
##           64 1 1 3 1 1 2 2 3 1  2  1  2  3  2  3  3  1  2  1  1  2
##           65 3 2 3 1 2 1 3 2 1  1  2  3  3  1  1  2  1  2  3  1  2
##           66 1 1 3 2 3 2 1 3 2  3  3  2  3  3  1  1  3  3  2  2  3
##           67 3 3 1 3 3 3 1 3 3  1  2  1  2  3  2  3  1  1  2  3  3
##           68 3 2 2 2 2 3 2 1 1  1  3  1  1  1  1  2  1  3  1  1  3
##           69 1 1 1 3 2 1 1 3 2  3  3  3  1  1  3  1  1  3  3  3  3
##           70 3 1 1 2 2 2 3 3 2  1  2  1  1  2  3  2  3  3  2  2  3
##           71 2 2 2 3 2 3 3 2 1  1  3  1  3  1  2  2  2  1  3  3  3
##           72 1 1 2 2 2 2 2 2 2  3  2  3  3  2  2  2  3  3  3  2  1
##           73 3 3 2 3 3 3 2 1 3  1  2  3  3  3  2  1  1  2  3  2  3
##           74 1 2 2 1 2 1 2 1 2  1  2  2  2  3  2  2  1  2  3  3  1
##           75 3 3 1 1 1 1 3 2 3  3  1  1  2  3  1  2  1  3  2  3  1
##           76 2 2 2 2 3 1 3 2 1  1  1  3  2  2  1  1  3  1  2  3  1
##           77 2 2 3 2 1 2 3 2 1  3  1  2  3  1  1  2  3  3  2  3  2
##           78 3 3 1 2 3 1 3 3 1  2  3  3  1  2  2  2  1  3  1  1  3
##           79 1 3 1 2 3 3 3 2 1  1  2  3  2  1  2  2  3  3  1  2  2
##           80 2 1 2 2 2 2 1 2 3  1  3  1  1  2  1  1  2  3  1  2  2
##           81 2 3 1 3 2 2 1 3 1  2  1  2  3  3  1  3  3  3  1  2  3
##           82 3 3 3 2 2 1 2 1 1  2  1  3  1  1  1  1  1  1  1  3  1
##           83 1 2 1 1 1 1 2 2 2  3  1  3  3  1  3  2  3  1  1  3  1
##           84 2 3 2 2 1 1 1 3 1  2  1  1  3  3  3  2  1  1  2  1  3
##           85 3 2 3 3 2 1 3 3 1  2  3  1  1  1  3  1  1  3  1  3  3
##           86 3 1 2 3 3 2 1 3 3  1  2  3  1  2  3  3  1  1  3  1  2
##           87 1 3 3 1 2 3 3 3 2  3  1  1  3  1  2  1  3  3  2  3  1
##           88 1 2 2 2 1 3 2 3 1  1  3  2  2  3  1  2  2  2  1  2  3
##           89 3 3 1 3 2 3 1 1 3  3  2  2  3  3  1  2  2  1  3  2  3
##           90 3 3 2 2 1 3 1 3 2  3  1  1  1  1  3  2  1  1  2  3  1
##           91 3 2 1 1 1 3 2 3 1  1  1  3  2  1  2  3  2  2  2  3  3
##           92 1 3 1 2 1 2 2 2 3  1  3  1  1  3  3  2  2  1  1  2  2
##           93 2 2 1 1 3 3 1 1 2  3  2  2  1  3  1  1  2  2  3  2  1
##           94 2 2 1 2 3 1 1 1 2  3  1  2  1  1  3  1  2  1  1  2  2
##           95 2 1 3 2 3 2 3 2 2  2  3  2  1  1  2  1  3  2  2  3  3
##           96 2 3 3 3 3 1 3 3 1  2  3  1  3  1  3  2  1  3  2  1  3
##           97 1 3 3 3 2 3 1 2 1  2  2  3  2  3  2  1  3  3  2  3  3
##           98 3 3 2 2 3 3 2 2 3  2  2  1  2  1  3  2  1  2  3  3  2
##           99 1 3 2 3 1 3 3 3 1  1  3  1  3  2  1  1  2  1  3  3  1
##          100 2 2 2 1 3 1 2 3 2  1  3  2  1  2  1  3  3  2  2  1  1
##          101 1 3 2 3 3 2 1 3 1  1  1  1  1  2  3  3  2  3  1  3  2
##          102 1 2 1 2 3 2 1 2 3  3  1  3  3  1  3  3  2  3  3  3  2
##          103 1 1 3 1 2 1 2 3 1  2  2  1  2  1  2  2  3  3  2  3  3
##          104 3 3 1 3 3 1 3 3 1  2  1  1  3  3  3  1  3  1  3  1  2
##          105 1 3 2 1 2 2 2 3 2  1  1  2  1  2  1  1  3  2  3  1  2
##          106 3 3 1 2 2 2 1 1 2  2  1  1  2  3  2  2  2  1  1  1  1
##          107 3 3 2 2 2 2 3 2 2  2  3  2  2  1  1  1  2  1  2  1  1
##          108 2 1 2 2 1 2 1 1 1  3  1  3  3  1  1  2  1  2  3  3  3
##          109 1 2 2 1 3 1 2 1 3  1  2  3  1  2  3  1  3  1  3  1  3
##          110 2 3 1 2 3 3 3 2 2  1  1  2  1  1  2  2  1  2  3  1  1
##          111 3 2 2 3 1 2 1 1 3  2  2  1  1  2  1  2  1  2  2  1  3
##          112 3 1 1 1 1 2 1 3 1  2  2  1  2  1  1  2  1  2  3  1  2
##          113 2 2 3 2 1 3 3 1 3  1  1  3  2  3  2  1  2  1  3  1  2
##          114 2 3 1 1 3 3 1 1 3  1  3  2  3  2  1  3  3  2  1  1  1
##          115 3 1 3 1 2 3 3 1 1  2  3  1  3  2  1  1  3  1  1  2  2
##          116 3 3 1 1 1 2 2 1 3  2  1  1  1  1  1  2  3  3  3  2  1
##          117 1 2 3 2 3 1 1 3 3  1  2  3  3  2  1  3  3  2  2  2  2
##          118 2 2 2 3 3 2 1 2 1  2  1  3  2  2  3  2  2  1  1  2  2
##          119 3 1 1 3 3 1 1 1 2  3  3  3  3  3  3  2  1  3  3  2  3
##          120 1 2 1 3 1 3 3 3 1  1  1  2  3  1  3  2  1  3  3  1  3

Check visual

# Sample data (replace this with your actual data)
# df <- data.frame(CommAtt = ..., Eff = ...)

# Create a boxplot comparing Effort across different Communicative Attempts
ggplot(final_data, aes(x = as.factor(CommAtt), y = Eff)) +
  geom_boxplot(aes(fill = as.factor(CommAtt))) +  # Boxplot with fill based on CommAtt
  labs(title = "Comparison of Effort Across Communicative Attempts",
       x = "Communicative Attempts",
       y = "Effort",
       fill = "CommAtt") + 
  theme_minimal() +
  theme(legend.position = "none")  # Optional: remove the legend

Modeling

Now let’s model our synth data

Contrast coding

# Convert necessary columns to factors
final_data$CommAtt <- as.factor(final_data$CommAtt)
final_data$Modality <- as.factor(final_data$Modality)
final_data$Participant <- as.factor(final_data$Participant)
final_data$Concept <- as.factor(final_data$Concept)
final_data$TrialNumber <- as.numeric(final_data$TrialNumber)  # Ensure TrialNumber is numeric

DP: Check contrasts of factors

contrasts(final_data$CommAtt) <- MASS::contr.sdif(3)
contrasts(final_data$Modality) <- contr.sum(3)/2

Center trial number

final_data$TrialNumber_c <- final_data$TrialNumber - median(range(final_data$TrialNumber))
range(final_data$TrialNumber_c)
## [1] -26.5  26.5
range(final_data$TrialNumber)
## [1]  1 54

Standardize (z-score) all continuous predictors.

But if the measures of familiarity, expressibility, and big5 are on a rating scale we can just subtract the median again from these to centre them (and we don’t need to standardize).

final_data$Familiarity <- final_data$Familiarity - median(range(final_data$Familiarity))
final_data$Big5 <- final_data$Big5 - median(range(final_data$Big5))

Z-score expressibility (because it’s continuous)

final_data <-
  final_data |> 
  mutate(Expressibility_z = (Expressibility - mean(Expressibility))/ sd(final_data$Expressibility, na.rm = T))

Simple linear mixed model

Before we do full Bayesian, we do quick LMM to see whether we are on good track

# Fit the linear mixed-effects model
lmer_model <- lmer(Eff ~ CommAtt + Familiarity + Big5 + Expressibility + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept), data = final_data)

# Summary of the model
summary(model)

# Check model diagnostics
# Plot residuals vs fitted values
plot(model)

# Extract coefficients
coefficients <- summary(model)$coefficients
print(coefficients)

# If you want to save the model output
saveRDS(model, here("06_Modeling", "models", "linear_mixed_effects_model.rds"))
lmer_model <- readRDS(here("06_Modeling", "models", "linear_mixed_effects_model.rds"))
summary(lmer_model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Eff ~ CommAtt + Familiarity + Big5 + Expressibility + TrialNumber +  
##     Modality + (1 | Participant) + (1 | Concept)
##    Data: final_data
## 
## REML criterion at convergence: 8420.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6533 -0.6248 -0.0245  0.6415  4.9419 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Participant (Intercept) 0.013366 0.11561 
##  Concept     (Intercept) 0.001346 0.03668 
##  Residual                0.297073 0.54504 
## Number of obs: 5067, groups:  Participant, 120; Concept, 21
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)      0.0748620  0.0449111   1.667
## CommAtt2         1.1706126  0.0171346  68.319
## CommAtt3         0.1163042  0.0217732   5.342
## Familiarity      1.1668413  0.0483836  24.116
## Big5             1.2090271  0.0460398  26.260
## Expressibility   1.1281766  0.0294578  38.298
## TrialNumber     -0.0011358  0.0006236  -1.821
## Modalitygesture  0.9456262  0.0203315  46.511
## Modalityvocal    0.9134778  0.0206897  44.151
## 
## Correlation of Fixed Effects:
##             (Intr) CmmAt2 CmmAt3 Fmlrty Big5   Exprss TrlNmb Mdltyg
## CommAtt2    -0.148                                                 
## CommAtt3    -0.104  0.320                                          
## Familiarity -0.567  0.003  0.001                                   
## Big5        -0.447  0.003 -0.001 -0.092                            
## Expressblty -0.387  0.004  0.006  0.014  0.004                     
## TrialNumber -0.288 -0.027 -0.044  0.006  0.005 -0.012              
## Modaltygstr -0.332 -0.010 -0.018  0.004  0.010  0.351 -0.008       
## Modalityvcl -0.325 -0.006 -0.014  0.005 -0.003  0.380 -0.047  0.562

The Betas seem quite close to how we create the synthetic data, probably there are some other moderations since the causal relations between the variables are quite complex

Open question: dyad + participant

How/whether to include participant in addition to dyad in the random effects

Our thoughts: participant could maybe be nested within dyad, or do we not need participant? Must be something suggested in the literature.

Participant probably must be included because the Big 5 personality trait measure corresponds to each participant within the dyad.

Check: https://github.com/FredericBlum/exploring-individual-variation-in-turkish-heritage-speakers-complex-linguistic-productions

Onur and Frede on nesting:

First, we do not use mono- versus bilingualism or even the group variable as a main independent variable in our model. Instead, it is incorporated in a nested random effect. This indicates that we do not view the speaker group status as a main driver of DM production. Second, the group variable that we utilize does not take one of the three speaker groups as a baseline. We sum-coded the variable which incorporates the idea that monolinguals cannot be an adequate baseline for bilinguals.

Check model diagnostics

performance::check_model(lmer_model)
## Failed to compute posterior predictive checks with `re_formula=NULL`.
##   Trying again with `re_formula=NA` now.

library(sjPlot)
## Warning: Paket 'sjPlot' wurde unter R Version 4.3.3 erstellt
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
sjPlot::plot_model(lmer_model)

library(ggplot2)
ggeffects::ggpredict(lmer_model,
                     terms = c("CommAtt"),
                     type = "random",
                     interval = "confidence") |> 
  as_tibble() |> 
  ggplot2::ggplot() +
  aes(x = x, y = predicted) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
  geom_line(group = 1) +
  theme_bw()

Going Bayesian

Potential ways:

However, after consultation with Daniela, we can get rid of the non-linear issue by contrast coding

library(ggplot2)
library(patchwork)
## Warning: Paket 'patchwork' wurde unter R Version 4.3.3 erstellt
library(bayesplot)
## Warning: Paket 'bayesplot' wurde unter R Version 4.3.3 erstellt
## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(brms)
## Warning: Paket 'brms' wurde unter R Version 4.3.3 erstellt
## Lade nötiges Paket: Rcpp
## Warning: Paket 'Rcpp' wurde unter R Version 4.3.3 erstellt
## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attache Paket: 'brms'
## Das folgende Objekt ist maskiert 'package:bayesplot':
## 
##     rhat
## Das folgende Objekt ist maskiert 'package:lme4':
## 
##     ngrps
## Das folgende Objekt ist maskiert 'package:stats':
## 
##     ar
library(beepr)
## Warning: Paket 'beepr' wurde unter R Version 4.3.3 erstellt
# use all available cores for parallel computing
options(mc.cores = parallel::detectCores())

Model 1: starter

Let’s start slowly, first with simple model accounting only for outcome + main predictor

Eff ~ Normal(μ,σ^2) μ= β0 + β1*CommAtt

fit_eff1 <- brm(Eff ~ 1 + CommAtt, data=final_data)
saveRDS(fit_eff1, here("06_Modeling", "models", "fit_eff1.rds"))

beep(5)
fit_eff1 <- readRDS(here("06_Modeling", "models", "fit_eff1.rds"))

summary(fit_eff1)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      2.74      0.01     2.72     2.77 1.00     4017     3144
## CommAtt2M1     1.58      0.03     1.52     1.63 1.00     3578     3162
## CommAtt3M2    -2.12      0.04    -2.19    -2.04 1.00     3399     2958
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.89      0.01     0.88     0.91 1.00     4138     3305
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# even so simple model can quite well capture the relationship in the simulated data

plot(conditional_effects(fit_eff1), points = TRUE)

# but there is quite a lot of variation, causing overlap

pp_check(fit_eff1, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# we see that the simulated data are not that loyal to the data

pp_check(fit_eff1, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# instead of cloud of points, we see almost perfectly linear relationship between outcome and residuals, indicating strong problems with the independence assumption of the errors.

Model 2: varying intercept (random effect)

Now we know we will probably need partial pooling for participants, because they share lot of commonalities (Later we might adapt this, because it might be participant within dyads. Dyads might also share lot of commonalities)

fit_eff2 <- brm(Eff ~ 1 + CommAtt + (1 | Participant), data=final_data)
saveRDS(fit_eff2, here("06_Modeling", "models", "fit_eff2.rds"))

beep(5)
fit_eff2 <- readRDS(here("06_Modeling", "models", "fit_eff2.rds"))

summary(fit_eff2)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + (1 | Participant) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Participant (Number of levels: 120) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.51      0.03     0.45     0.58 1.01      406      795
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      2.75      0.05     2.65     2.84 1.02      204      383
## CommAtt2M1     1.59      0.02     1.54     1.63 1.00     4096     3460
## CommAtt3M2    -2.12      0.03    -2.18    -2.06 1.00     3988     3397
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.74      0.01     0.72     0.75 1.00     4894     3154
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(conditional_effects(fit_eff2), points = TRUE)

# stays the same

pp_check(fit_eff2, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# not really improvement

pp_check(fit_eff2, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# ok, here we see some improvement here!

# checking to coeff of the parameters
coef(fit_eff2)
## $Participant
## , , Intercept
## 
##     Estimate  Est.Error     Q2.5    Q97.5
## 1   3.579610 0.11441565 3.352994 3.812583
## 2   2.522447 0.11044733 2.314125 2.741480
## 3   3.146291 0.10402491 2.944224 3.347758
## 4   2.262994 0.10891029 2.051830 2.477337
## 5   1.668679 0.11156294 1.453740 1.885679
## 6   3.210786 0.10755374 3.001736 3.421982
## 7   1.743726 0.09956315 1.546578 1.942408
## 8   2.415176 0.11570030 2.188098 2.644443
## 9   2.941818 0.10847545 2.726732 3.153139
## 10  2.850522 0.11247866 2.631404 3.071326
## 11  3.235421 0.10197221 3.035279 3.431147
## 12  3.268229 0.10853672 3.055044 3.486600
## 13  3.101527 0.11377116 2.872792 3.323975
## 14  2.134733 0.10677733 1.927514 2.340458
## 15  2.820435 0.11210555 2.603776 3.041106
## 16  2.327921 0.10479141 2.120863 2.528586
## 17  2.612322 0.10537909 2.405861 2.819664
## 18  3.023965 0.11462556 2.795903 3.246746
## 19  2.778095 0.11530711 2.551082 3.002631
## 20  3.151138 0.10798837 2.942186 3.357480
## 21  3.567949 0.11302305 3.352112 3.786573
## 22  2.429940 0.11628572 2.196410 2.659209
## 23  2.968761 0.11250622 2.749033 3.188065
## 24  2.704271 0.10887674 2.491699 2.912234
## 25  2.672293 0.10259883 2.471873 2.872144
## 26  2.779301 0.10857503 2.573211 2.993449
## 27  2.465089 0.11036437 2.249112 2.681764
## 28  2.223275 0.11764174 1.989483 2.456018
## 29  1.846008 0.11587214 1.628159 2.072764
## 30  3.048777 0.11053537 2.831188 3.261401
## 31  2.809188 0.10616442 2.601615 3.014372
## 32  2.543952 0.10470197 2.338840 2.751140
## 33  2.552387 0.11051521 2.331680 2.763790
## 34  2.076140 0.11182684 1.856244 2.292397
## 35  2.219996 0.11277154 1.995945 2.436680
## 36  2.493701 0.10603293 2.281718 2.702046
## 37  3.679949 0.11266136 3.461175 3.900210
## 38  2.736239 0.11295006 2.512656 2.950943
## 39  2.702837 0.11809388 2.469696 2.933019
## 40  2.882159 0.10722781 2.678426 3.093300
## 41  2.472661 0.11187430 2.253757 2.687813
## 42  2.745807 0.11707738 2.510145 2.974322
## 43  2.702553 0.10588311 2.497378 2.912475
## 44  2.260898 0.11366341 2.040076 2.493751
## 45  3.384675 0.11866167 3.152057 3.616854
## 46  2.716896 0.10449897 2.511432 2.925632
## 47  1.852768 0.10364244 1.652043 2.051654
## 48  3.342887 0.11549076 3.118391 3.559941
## 49  1.771033 0.11031665 1.555702 1.992023
## 50  3.346107 0.12714312 3.099476 3.597375
## 51  3.877450 0.11178799 3.661954 4.090479
## 52  1.829994 0.10348166 1.624729 2.030700
## 53  2.396591 0.11912672 2.162229 2.638269
## 54  2.380222 0.11579929 2.150989 2.602008
## 55  3.182091 0.10873405 2.974668 3.395113
## 56  2.315297 0.11311378 2.096352 2.537675
## 57  2.064963 0.11308449 1.837772 2.284046
## 58  2.314193 0.11107981 2.098013 2.523934
## 59  2.399026 0.11424194 2.180342 2.626832
## 60  2.089906 0.11778003 1.866431 2.316164
## 61  2.700226 0.11952011 2.469309 2.926424
## 62  2.619516 0.12051327 2.384422 2.857954
## 63  3.171361 0.11073947 2.951481 3.382997
## 64  2.611329 0.11696504 2.383084 2.841700
## 65  2.688010 0.11691878 2.461584 2.922726
## 66  3.127317 0.10673919 2.919742 3.336950
## 67  2.255830 0.10563610 2.052473 2.461304
## 68  2.857681 0.12025854 2.624437 3.091862
## 69  3.289623 0.11099989 3.076365 3.506857
## 70  2.945625 0.10547291 2.740275 3.152410
## 71  3.512916 0.10935253 3.291953 3.730189
## 72  2.990230 0.10685282 2.780816 3.191803
## 73  1.771854 0.10382751 1.571826 1.971438
## 74  2.756099 0.11608226 2.531082 2.981185
## 75  3.400385 0.11514800 3.176591 3.622344
## 76  2.866797 0.11839811 2.632107 3.091283
## 77  2.640661 0.11011714 2.426721 2.854269
## 78  2.335450 0.10309816 2.133630 2.532908
## 79  3.561743 0.10861026 3.351745 3.772179
## 80  3.247257 0.11588270 3.016624 3.480770
## 81  2.780556 0.10592433 2.576949 2.991000
## 82  1.933096 0.12234709 1.691496 2.181393
## 83  2.621529 0.11831408 2.391282 2.850984
## 84  3.181666 0.11442281 2.960273 3.404832
## 85  3.313460 0.10981852 3.101184 3.523723
## 86  3.107441 0.10816763 2.899903 3.327315
## 87  2.968963 0.10515426 2.768879 3.171367
## 88  2.946239 0.11383147 2.728444 3.167263
## 89  2.801226 0.10178248 2.600620 2.999539
## 90  3.110846 0.11296329 2.897673 3.326217
## 91  2.349701 0.11047740 2.129544 2.561816
## 92  3.065909 0.11395658 2.840459 3.290479
## 93  3.123400 0.11582261 2.903462 3.350258
## 94  2.696007 0.12183158 2.459617 2.932619
## 95  2.789953 0.10861537 2.574905 3.002887
## 96  3.029803 0.10494297 2.822342 3.237188
## 97  3.239874 0.10240338 3.044348 3.440797
## 98  2.175207 0.10657819 1.972430 2.380844
## 99  3.136929 0.11506284 2.914134 3.358832
## 100 1.889527 0.10989286 1.671083 2.102376
## 101 3.043525 0.11398017 2.819380 3.263130
## 102 2.275650 0.10904142 2.071340 2.491571
## 103 2.072891 0.10972387 1.856593 2.286583
## 104 2.762111 0.10780242 2.551544 2.978329
## 105 1.945242 0.11819705 1.708879 2.172253
## 106 3.111959 0.12023253 2.879270 3.347590
## 107 2.747825 0.11557633 2.528067 2.970314
## 108 2.987815 0.11518119 2.764853 3.213802
## 109 3.613064 0.11517567 3.383850 3.836056
## 110 3.726170 0.11354328 3.502244 3.946415
## 111 2.775667 0.11775600 2.542527 3.008275
## 112 2.727688 0.12697156 2.477013 2.969186
## 113 3.086444 0.10924760 2.868714 3.301903
## 114 3.506301 0.11439792 3.281420 3.723125
## 115 2.603470 0.11288770 2.385401 2.817476
## 116 2.761380 0.11791953 2.526261 2.993304
## 117 3.276696 0.10850088 3.067212 3.495088
## 118 2.783191 0.11433916 2.558433 3.005766
## 119 3.048421 0.10547291 2.835668 3.258379
## 120 1.960899 0.10806763 1.747192 2.173225
## 
## , , CommAtt2M1
## 
##     Estimate  Est.Error     Q2.5    Q97.5
## 1   1.588552 0.02323773 1.542844 1.633812
## 2   1.588552 0.02323773 1.542844 1.633812
## 3   1.588552 0.02323773 1.542844 1.633812
## 4   1.588552 0.02323773 1.542844 1.633812
## 5   1.588552 0.02323773 1.542844 1.633812
## 6   1.588552 0.02323773 1.542844 1.633812
## 7   1.588552 0.02323773 1.542844 1.633812
## 8   1.588552 0.02323773 1.542844 1.633812
## 9   1.588552 0.02323773 1.542844 1.633812
## 10  1.588552 0.02323773 1.542844 1.633812
## 11  1.588552 0.02323773 1.542844 1.633812
## 12  1.588552 0.02323773 1.542844 1.633812
## 13  1.588552 0.02323773 1.542844 1.633812
## 14  1.588552 0.02323773 1.542844 1.633812
## 15  1.588552 0.02323773 1.542844 1.633812
## 16  1.588552 0.02323773 1.542844 1.633812
## 17  1.588552 0.02323773 1.542844 1.633812
## 18  1.588552 0.02323773 1.542844 1.633812
## 19  1.588552 0.02323773 1.542844 1.633812
## 20  1.588552 0.02323773 1.542844 1.633812
## 21  1.588552 0.02323773 1.542844 1.633812
## 22  1.588552 0.02323773 1.542844 1.633812
## 23  1.588552 0.02323773 1.542844 1.633812
## 24  1.588552 0.02323773 1.542844 1.633812
## 25  1.588552 0.02323773 1.542844 1.633812
## 26  1.588552 0.02323773 1.542844 1.633812
## 27  1.588552 0.02323773 1.542844 1.633812
## 28  1.588552 0.02323773 1.542844 1.633812
## 29  1.588552 0.02323773 1.542844 1.633812
## 30  1.588552 0.02323773 1.542844 1.633812
## 31  1.588552 0.02323773 1.542844 1.633812
## 32  1.588552 0.02323773 1.542844 1.633812
## 33  1.588552 0.02323773 1.542844 1.633812
## 34  1.588552 0.02323773 1.542844 1.633812
## 35  1.588552 0.02323773 1.542844 1.633812
## 36  1.588552 0.02323773 1.542844 1.633812
## 37  1.588552 0.02323773 1.542844 1.633812
## 38  1.588552 0.02323773 1.542844 1.633812
## 39  1.588552 0.02323773 1.542844 1.633812
## 40  1.588552 0.02323773 1.542844 1.633812
## 41  1.588552 0.02323773 1.542844 1.633812
## 42  1.588552 0.02323773 1.542844 1.633812
## 43  1.588552 0.02323773 1.542844 1.633812
## 44  1.588552 0.02323773 1.542844 1.633812
## 45  1.588552 0.02323773 1.542844 1.633812
## 46  1.588552 0.02323773 1.542844 1.633812
## 47  1.588552 0.02323773 1.542844 1.633812
## 48  1.588552 0.02323773 1.542844 1.633812
## 49  1.588552 0.02323773 1.542844 1.633812
## 50  1.588552 0.02323773 1.542844 1.633812
## 51  1.588552 0.02323773 1.542844 1.633812
## 52  1.588552 0.02323773 1.542844 1.633812
## 53  1.588552 0.02323773 1.542844 1.633812
## 54  1.588552 0.02323773 1.542844 1.633812
## 55  1.588552 0.02323773 1.542844 1.633812
## 56  1.588552 0.02323773 1.542844 1.633812
## 57  1.588552 0.02323773 1.542844 1.633812
## 58  1.588552 0.02323773 1.542844 1.633812
## 59  1.588552 0.02323773 1.542844 1.633812
## 60  1.588552 0.02323773 1.542844 1.633812
## 61  1.588552 0.02323773 1.542844 1.633812
## 62  1.588552 0.02323773 1.542844 1.633812
## 63  1.588552 0.02323773 1.542844 1.633812
## 64  1.588552 0.02323773 1.542844 1.633812
## 65  1.588552 0.02323773 1.542844 1.633812
## 66  1.588552 0.02323773 1.542844 1.633812
## 67  1.588552 0.02323773 1.542844 1.633812
## 68  1.588552 0.02323773 1.542844 1.633812
## 69  1.588552 0.02323773 1.542844 1.633812
## 70  1.588552 0.02323773 1.542844 1.633812
## 71  1.588552 0.02323773 1.542844 1.633812
## 72  1.588552 0.02323773 1.542844 1.633812
## 73  1.588552 0.02323773 1.542844 1.633812
## 74  1.588552 0.02323773 1.542844 1.633812
## 75  1.588552 0.02323773 1.542844 1.633812
## 76  1.588552 0.02323773 1.542844 1.633812
## 77  1.588552 0.02323773 1.542844 1.633812
## 78  1.588552 0.02323773 1.542844 1.633812
## 79  1.588552 0.02323773 1.542844 1.633812
## 80  1.588552 0.02323773 1.542844 1.633812
## 81  1.588552 0.02323773 1.542844 1.633812
## 82  1.588552 0.02323773 1.542844 1.633812
## 83  1.588552 0.02323773 1.542844 1.633812
## 84  1.588552 0.02323773 1.542844 1.633812
## 85  1.588552 0.02323773 1.542844 1.633812
## 86  1.588552 0.02323773 1.542844 1.633812
## 87  1.588552 0.02323773 1.542844 1.633812
## 88  1.588552 0.02323773 1.542844 1.633812
## 89  1.588552 0.02323773 1.542844 1.633812
## 90  1.588552 0.02323773 1.542844 1.633812
## 91  1.588552 0.02323773 1.542844 1.633812
## 92  1.588552 0.02323773 1.542844 1.633812
## 93  1.588552 0.02323773 1.542844 1.633812
## 94  1.588552 0.02323773 1.542844 1.633812
## 95  1.588552 0.02323773 1.542844 1.633812
## 96  1.588552 0.02323773 1.542844 1.633812
## 97  1.588552 0.02323773 1.542844 1.633812
## 98  1.588552 0.02323773 1.542844 1.633812
## 99  1.588552 0.02323773 1.542844 1.633812
## 100 1.588552 0.02323773 1.542844 1.633812
## 101 1.588552 0.02323773 1.542844 1.633812
## 102 1.588552 0.02323773 1.542844 1.633812
## 103 1.588552 0.02323773 1.542844 1.633812
## 104 1.588552 0.02323773 1.542844 1.633812
## 105 1.588552 0.02323773 1.542844 1.633812
## 106 1.588552 0.02323773 1.542844 1.633812
## 107 1.588552 0.02323773 1.542844 1.633812
## 108 1.588552 0.02323773 1.542844 1.633812
## 109 1.588552 0.02323773 1.542844 1.633812
## 110 1.588552 0.02323773 1.542844 1.633812
## 111 1.588552 0.02323773 1.542844 1.633812
## 112 1.588552 0.02323773 1.542844 1.633812
## 113 1.588552 0.02323773 1.542844 1.633812
## 114 1.588552 0.02323773 1.542844 1.633812
## 115 1.588552 0.02323773 1.542844 1.633812
## 116 1.588552 0.02323773 1.542844 1.633812
## 117 1.588552 0.02323773 1.542844 1.633812
## 118 1.588552 0.02323773 1.542844 1.633812
## 119 1.588552 0.02323773 1.542844 1.633812
## 120 1.588552 0.02323773 1.542844 1.633812
## 
## , , CommAtt3M2
## 
##      Estimate  Est.Error      Q2.5     Q97.5
## 1   -2.120611 0.03038427 -2.180229 -2.061068
## 2   -2.120611 0.03038427 -2.180229 -2.061068
## 3   -2.120611 0.03038427 -2.180229 -2.061068
## 4   -2.120611 0.03038427 -2.180229 -2.061068
## 5   -2.120611 0.03038427 -2.180229 -2.061068
## 6   -2.120611 0.03038427 -2.180229 -2.061068
## 7   -2.120611 0.03038427 -2.180229 -2.061068
## 8   -2.120611 0.03038427 -2.180229 -2.061068
## 9   -2.120611 0.03038427 -2.180229 -2.061068
## 10  -2.120611 0.03038427 -2.180229 -2.061068
## 11  -2.120611 0.03038427 -2.180229 -2.061068
## 12  -2.120611 0.03038427 -2.180229 -2.061068
## 13  -2.120611 0.03038427 -2.180229 -2.061068
## 14  -2.120611 0.03038427 -2.180229 -2.061068
## 15  -2.120611 0.03038427 -2.180229 -2.061068
## 16  -2.120611 0.03038427 -2.180229 -2.061068
## 17  -2.120611 0.03038427 -2.180229 -2.061068
## 18  -2.120611 0.03038427 -2.180229 -2.061068
## 19  -2.120611 0.03038427 -2.180229 -2.061068
## 20  -2.120611 0.03038427 -2.180229 -2.061068
## 21  -2.120611 0.03038427 -2.180229 -2.061068
## 22  -2.120611 0.03038427 -2.180229 -2.061068
## 23  -2.120611 0.03038427 -2.180229 -2.061068
## 24  -2.120611 0.03038427 -2.180229 -2.061068
## 25  -2.120611 0.03038427 -2.180229 -2.061068
## 26  -2.120611 0.03038427 -2.180229 -2.061068
## 27  -2.120611 0.03038427 -2.180229 -2.061068
## 28  -2.120611 0.03038427 -2.180229 -2.061068
## 29  -2.120611 0.03038427 -2.180229 -2.061068
## 30  -2.120611 0.03038427 -2.180229 -2.061068
## 31  -2.120611 0.03038427 -2.180229 -2.061068
## 32  -2.120611 0.03038427 -2.180229 -2.061068
## 33  -2.120611 0.03038427 -2.180229 -2.061068
## 34  -2.120611 0.03038427 -2.180229 -2.061068
## 35  -2.120611 0.03038427 -2.180229 -2.061068
## 36  -2.120611 0.03038427 -2.180229 -2.061068
## 37  -2.120611 0.03038427 -2.180229 -2.061068
## 38  -2.120611 0.03038427 -2.180229 -2.061068
## 39  -2.120611 0.03038427 -2.180229 -2.061068
## 40  -2.120611 0.03038427 -2.180229 -2.061068
## 41  -2.120611 0.03038427 -2.180229 -2.061068
## 42  -2.120611 0.03038427 -2.180229 -2.061068
## 43  -2.120611 0.03038427 -2.180229 -2.061068
## 44  -2.120611 0.03038427 -2.180229 -2.061068
## 45  -2.120611 0.03038427 -2.180229 -2.061068
## 46  -2.120611 0.03038427 -2.180229 -2.061068
## 47  -2.120611 0.03038427 -2.180229 -2.061068
## 48  -2.120611 0.03038427 -2.180229 -2.061068
## 49  -2.120611 0.03038427 -2.180229 -2.061068
## 50  -2.120611 0.03038427 -2.180229 -2.061068
## 51  -2.120611 0.03038427 -2.180229 -2.061068
## 52  -2.120611 0.03038427 -2.180229 -2.061068
## 53  -2.120611 0.03038427 -2.180229 -2.061068
## 54  -2.120611 0.03038427 -2.180229 -2.061068
## 55  -2.120611 0.03038427 -2.180229 -2.061068
## 56  -2.120611 0.03038427 -2.180229 -2.061068
## 57  -2.120611 0.03038427 -2.180229 -2.061068
## 58  -2.120611 0.03038427 -2.180229 -2.061068
## 59  -2.120611 0.03038427 -2.180229 -2.061068
## 60  -2.120611 0.03038427 -2.180229 -2.061068
## 61  -2.120611 0.03038427 -2.180229 -2.061068
## 62  -2.120611 0.03038427 -2.180229 -2.061068
## 63  -2.120611 0.03038427 -2.180229 -2.061068
## 64  -2.120611 0.03038427 -2.180229 -2.061068
## 65  -2.120611 0.03038427 -2.180229 -2.061068
## 66  -2.120611 0.03038427 -2.180229 -2.061068
## 67  -2.120611 0.03038427 -2.180229 -2.061068
## 68  -2.120611 0.03038427 -2.180229 -2.061068
## 69  -2.120611 0.03038427 -2.180229 -2.061068
## 70  -2.120611 0.03038427 -2.180229 -2.061068
## 71  -2.120611 0.03038427 -2.180229 -2.061068
## 72  -2.120611 0.03038427 -2.180229 -2.061068
## 73  -2.120611 0.03038427 -2.180229 -2.061068
## 74  -2.120611 0.03038427 -2.180229 -2.061068
## 75  -2.120611 0.03038427 -2.180229 -2.061068
## 76  -2.120611 0.03038427 -2.180229 -2.061068
## 77  -2.120611 0.03038427 -2.180229 -2.061068
## 78  -2.120611 0.03038427 -2.180229 -2.061068
## 79  -2.120611 0.03038427 -2.180229 -2.061068
## 80  -2.120611 0.03038427 -2.180229 -2.061068
## 81  -2.120611 0.03038427 -2.180229 -2.061068
## 82  -2.120611 0.03038427 -2.180229 -2.061068
## 83  -2.120611 0.03038427 -2.180229 -2.061068
## 84  -2.120611 0.03038427 -2.180229 -2.061068
## 85  -2.120611 0.03038427 -2.180229 -2.061068
## 86  -2.120611 0.03038427 -2.180229 -2.061068
## 87  -2.120611 0.03038427 -2.180229 -2.061068
## 88  -2.120611 0.03038427 -2.180229 -2.061068
## 89  -2.120611 0.03038427 -2.180229 -2.061068
## 90  -2.120611 0.03038427 -2.180229 -2.061068
## 91  -2.120611 0.03038427 -2.180229 -2.061068
## 92  -2.120611 0.03038427 -2.180229 -2.061068
## 93  -2.120611 0.03038427 -2.180229 -2.061068
## 94  -2.120611 0.03038427 -2.180229 -2.061068
## 95  -2.120611 0.03038427 -2.180229 -2.061068
## 96  -2.120611 0.03038427 -2.180229 -2.061068
## 97  -2.120611 0.03038427 -2.180229 -2.061068
## 98  -2.120611 0.03038427 -2.180229 -2.061068
## 99  -2.120611 0.03038427 -2.180229 -2.061068
## 100 -2.120611 0.03038427 -2.180229 -2.061068
## 101 -2.120611 0.03038427 -2.180229 -2.061068
## 102 -2.120611 0.03038427 -2.180229 -2.061068
## 103 -2.120611 0.03038427 -2.180229 -2.061068
## 104 -2.120611 0.03038427 -2.180229 -2.061068
## 105 -2.120611 0.03038427 -2.180229 -2.061068
## 106 -2.120611 0.03038427 -2.180229 -2.061068
## 107 -2.120611 0.03038427 -2.180229 -2.061068
## 108 -2.120611 0.03038427 -2.180229 -2.061068
## 109 -2.120611 0.03038427 -2.180229 -2.061068
## 110 -2.120611 0.03038427 -2.180229 -2.061068
## 111 -2.120611 0.03038427 -2.180229 -2.061068
## 112 -2.120611 0.03038427 -2.180229 -2.061068
## 113 -2.120611 0.03038427 -2.180229 -2.061068
## 114 -2.120611 0.03038427 -2.180229 -2.061068
## 115 -2.120611 0.03038427 -2.180229 -2.061068
## 116 -2.120611 0.03038427 -2.180229 -2.061068
## 117 -2.120611 0.03038427 -2.180229 -2.061068
## 118 -2.120611 0.03038427 -2.180229 -2.061068
## 119 -2.120611 0.03038427 -2.180229 -2.061068
## 120 -2.120611 0.03038427 -2.180229 -2.061068

Some more diagnostics

vars <- c("b_Intercept", "r_Participant[1,Intercept]", "r_Participant[2,Intercept]")
pairs(fit_eff2, variable = vars)

We see negative correlation between b0 and b0j (and b0j have weak positive correlation with one another)

draws_eff2 <- as.matrix(fit_eff2, variable = vars)
colnames(draws_eff2) <- c("b0", "b0_1", "b0_2")
cors_eff2 <- cor(draws_eff2)
print(cors_eff2, digits = 2)
##         b0  b0_1  b0_2
## b0    1.00 -0.36 -0.39
## b0_1 -0.36  1.00  0.12
## b0_2 -0.39  0.12  1.00

Correlation between the parameters

coef_eff2 <- coef(fit_eff2, summary = FALSE)

print(cor(coef_eff2$Participant[, 1:5, "Intercept"]), digits = 2)
##        1       2      3     4       5
## 1  1.000 -0.0178 0.0461 0.012 -0.0100
## 2 -0.018  1.0000 0.0039 0.063 -0.0084
## 3  0.046  0.0039 1.0000 0.057  0.0135
## 4  0.012  0.0626 0.0571 1.000  0.0378
## 5 -0.010 -0.0084 0.0135 0.038  1.0000

Per-subject predicitons into the plot

conditions <- make_conditions(final_data, "Participant")
ce <- conditional_effects(
fit_eff2, conditions = conditions,
re_formula = NULL
)
plot(ce, ncol = 6, points = TRUE)

# ok this is too many pcns to see anything:)

Compare the first two models

loo_eff1 <- loo(fit_eff1)
loo_eff2 <- loo(fit_eff2)

loo_compare(loo_eff1, loo_eff2)
##          elpd_diff se_diff
## fit_eff2    0.0       0.0 
## fit_eff1 -912.1      38.1

Model 3: varying slope

Let’s also add varying slope, because ComAtt will too vary across subjects

fit_eff3 <- brm(Eff ~ 1 + CommAtt + (1 + CommAtt | Participant), # note that there was a mistake before
                data=final_data,
                # Number of cores to use
                cores = 4,)
saveRDS(fit_eff3, here("06_Modeling", "models", "fit_eff3.rds"))

beep(5)
fit_eff3 <- readRDS(here("06_Modeling", "models", "fit_eff3.rds"))

summary(fit_eff3)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 1278 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + (1 | CommAtt + Participant) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~CommAtt (Number of levels: 3) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.72      1.94     0.08     5.71 1.27       11       55
## 
## ~Participant (Number of levels: 120) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.53      0.03     0.45     0.59 1.16       19      436
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      2.23      1.69     0.05     6.82 1.50        7       32
## CommAtt2M1     2.10      2.46    -3.39     6.78 1.15       61      668
## CommAtt3M2    -1.25      2.86    -7.81     3.37 1.14       20      106
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.74      0.01     0.72     0.75 1.41      557     1789
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(conditional_effects(fit_eff3), points = TRUE)

# now the condidence intervals are much wider

pp_check(fit_eff3, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# not really improvement

pp_check(fit_eff3, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# blobby, but still very correlated

Model 4: full DAG

Now let’s finally replicate our DAG

fit_eff4 <- brm(Eff ~ 1 + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept),
                data = final_data,
                cores = 4)

saveRDS(fit_eff4, here("06_Modeling", "models", "fit_eff4.rds"))

beep(5)
fit_eff4 <- readRDS(here("06_Modeling", "models", "fit_eff4.rds"))

# summary
summary(fit_eff4)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.03      0.01     0.00     0.06 1.00     1309     1353
## 
## ~Participant (Number of levels: 120) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.02      0.01     0.00     0.04 1.00     1886     2247
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            2.71      0.01     2.69     2.74 1.00     4845     3442
## CommAtt2M1           1.60      0.02     1.56     1.64 1.00     6017     3209
## CommAtt3M2          -2.12      0.03    -2.17    -2.07 1.00     6121     3464
## Familiarity          1.23      0.03     1.16     1.29 1.00     5385     3222
## Big5                 1.19      0.03     1.13     1.25 1.00     5173     3189
## Expressibility_z     0.39      0.01     0.37     0.41 1.00     4850     2713
## TrialNumber_c       -0.00      0.00    -0.00     0.00 1.00     3868     2355
## Modality1           -0.55      0.03    -0.60    -0.49 1.00     5645     2824
## Modality2            0.30      0.03     0.25     0.35 1.00     5678     3188
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.65      0.01     0.63     0.66 1.00     6786     2916
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# the coefficients again look similar to what we set in the data (except perhaps expressibility?)

plot(fit_eff4)

plot(conditional_effects(fit_eff4), points = TRUE)

# the CrI are now again in sensible width

pp_check(fit_eff4, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# not really improvement

pp_check(fit_eff4, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# blobby, but still correlated

# note that we are still missing proper non-def priors

Model 5: full + varying intercepts

fit_eff5 <- brm(Eff ~ 1 + CommAtt + Modality + (1 | Participant) + (1 | Concept) + (1 | Familiarity) + (1 | Big5) + (1 | Expressibility_z) + (1 | TrialNumber_c), 
                data = final_data,
                cores = 4)

# NOTE: I accidentally rewrote fit5 with model fit6 so this will need to be re-run and re-saved

saveRDS(fit_eff5, here("06_Modeling", "models", "fit_eff5.rds"))

beep(5)
fit_eff5 <- readRDS(here("06_Modeling", "models", "fit_eff5.rds"))

# summary
summary(fit_eff5)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 41 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.06      0.05     0.00     0.17 1.02      140
## sd(CommAtt2M1)                 0.03      0.03     0.00     0.10 1.02      403
## sd(CommAtt3M2)                 0.03      0.03     0.00     0.11 1.01      361
## cor(Intercept,CommAtt2M1)      0.03      0.51    -0.87     0.91 1.01      874
## cor(Intercept,CommAtt3M2)     -0.15      0.52    -0.94     0.83 1.01      616
## cor(CommAtt2M1,CommAtt3M2)    -0.18      0.52    -0.95     0.84 1.01      632
##                            Tail_ESS
## sd(Intercept)                   104
## sd(CommAtt2M1)                  151
## sd(CommAtt3M2)                  131
## cor(Intercept,CommAtt2M1)      1145
## cor(Intercept,CommAtt3M2)       443
## cor(CommAtt2M1,CommAtt3M2)      286
## 
## ~Concept (Number of levels: 21) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.06      0.04     0.00     0.13 1.01      338     1047
## 
## ~Expressibility_z (Number of levels: 63) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.49      0.07     0.38     0.63 1.03      170
## sd(CommAtt2M1)                 0.61      0.06     0.51     0.72 1.01      624
## sd(CommAtt3M2)                 0.35      0.04     0.28     0.43 1.01      779
## cor(Intercept,CommAtt2M1)      0.92      0.06     0.79     0.98 1.02      128
## cor(Intercept,CommAtt3M2)     -0.95      0.04    -1.00    -0.86 1.01      294
## cor(CommAtt2M1,CommAtt3M2)    -0.81      0.06    -0.91    -0.68 1.00     1107
##                            Tail_ESS
## sd(Intercept)                   261
## sd(CommAtt2M1)                 1230
## sd(CommAtt3M2)                 1843
## cor(Intercept,CommAtt2M1)       172
## cor(Intercept,CommAtt3M2)       614
## cor(CommAtt2M1,CommAtt3M2)     1761
## 
## ~Familiarity (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.16      0.18     0.00     0.51 1.56        7
## sd(CommAtt2M1)                 0.09      0.11     0.00     0.32 1.53        7
## sd(CommAtt3M2)                 0.12      0.15     0.00     0.43 1.53        7
## cor(Intercept,CommAtt2M1)      0.26      0.59    -0.83     0.99 1.50        7
## cor(Intercept,CommAtt3M2)     -0.38      0.55    -1.00     0.78 1.51        7
## cor(CommAtt2M1,CommAtt3M2)    -0.40      0.55    -1.00     0.77 1.51        7
##                            Tail_ESS
## sd(Intercept)                    39
## sd(CommAtt2M1)                   29
## sd(CommAtt3M2)                   29
## cor(Intercept,CommAtt2M1)        39
## cor(Intercept,CommAtt3M2)        39
## cor(CommAtt2M1,CommAtt3M2)       32
## 
## ~Participant (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.36      0.18     0.01     0.53 1.54        7
## sd(CommAtt2M1)                 0.23      0.12     0.00     0.34 1.53        7
## sd(CommAtt3M2)                 0.30      0.16     0.00     0.45 1.53        7
## cor(Intercept,CommAtt2M1)      0.71      0.49    -0.70     0.99 1.51        7
## cor(Intercept,CommAtt3M2)     -0.76      0.44    -1.00     0.57 1.52        7
## cor(CommAtt2M1,CommAtt3M2)    -0.78      0.43    -1.00     0.57 1.52        7
##                            Tail_ESS
## sd(Intercept)                    27
## sd(CommAtt2M1)                   31
## sd(CommAtt3M2)                   29
## cor(Intercept,CommAtt2M1)        29
## cor(Intercept,CommAtt3M2)        28
## cor(CommAtt2M1,CommAtt3M2)       28
## 
## ~TrialNumber_c (Number of levels: 54) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.02      0.01     0.00     0.05 1.00     1085     1283
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      2.76      0.07     2.61     2.90 1.00      361      650
## CommAtt2M1     1.62      0.08     1.46     1.78 1.00      414      698
## CommAtt3M2    -2.15      0.06    -2.26    -2.02 1.01      417      808
## Modality1     -0.87      0.20    -1.24    -0.49 1.02      129      123
## Modality2      0.45      0.11     0.23     0.66 1.02      190      235
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.56      0.01     0.54     0.57 1.00     3809     1448
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# the coefficients again look similar to what we set in the data (except perhaps expressibility?)

plot(fit_eff5)

plot(conditional_effects(fit_eff5), points = TRUE)

# the CrI are now again in sensible width

pp_check(fit_eff5, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# not really improvement

pp_check(fit_eff5, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# blobby, but still correlated

Model 6: full with varying slopes and intercepts

So we not only assume that different participants have different base level for effort. But also that they differently ‘response’ to the communicative attempt

Same for familiarity, Big5, expressibility

fit_eff6 <- brm(Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c), 
                data = final_data,
                cores = 4)

saveRDS(fit_eff6, here("06_Modeling", "models", "fit_eff6.rds"))
beep(5)
fit_eff6 <- readRDS(here("06_Modeling", "models", "fit_eff6.rds"))

# summary
summary(fit_eff6)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 41 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.06      0.05     0.00     0.17 1.02      140
## sd(CommAtt2M1)                 0.03      0.03     0.00     0.10 1.02      403
## sd(CommAtt3M2)                 0.03      0.03     0.00     0.11 1.01      361
## cor(Intercept,CommAtt2M1)      0.03      0.51    -0.87     0.91 1.01      874
## cor(Intercept,CommAtt3M2)     -0.15      0.52    -0.94     0.83 1.01      616
## cor(CommAtt2M1,CommAtt3M2)    -0.18      0.52    -0.95     0.84 1.01      632
##                            Tail_ESS
## sd(Intercept)                   104
## sd(CommAtt2M1)                  151
## sd(CommAtt3M2)                  131
## cor(Intercept,CommAtt2M1)      1145
## cor(Intercept,CommAtt3M2)       443
## cor(CommAtt2M1,CommAtt3M2)      286
## 
## ~Concept (Number of levels: 21) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.06      0.04     0.00     0.13 1.01      338     1047
## 
## ~Expressibility_z (Number of levels: 63) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.49      0.07     0.38     0.63 1.03      170
## sd(CommAtt2M1)                 0.61      0.06     0.51     0.72 1.01      624
## sd(CommAtt3M2)                 0.35      0.04     0.28     0.43 1.01      779
## cor(Intercept,CommAtt2M1)      0.92      0.06     0.79     0.98 1.02      128
## cor(Intercept,CommAtt3M2)     -0.95      0.04    -1.00    -0.86 1.01      294
## cor(CommAtt2M1,CommAtt3M2)    -0.81      0.06    -0.91    -0.68 1.00     1107
##                            Tail_ESS
## sd(Intercept)                   261
## sd(CommAtt2M1)                 1230
## sd(CommAtt3M2)                 1843
## cor(Intercept,CommAtt2M1)       172
## cor(Intercept,CommAtt3M2)       614
## cor(CommAtt2M1,CommAtt3M2)     1761
## 
## ~Familiarity (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.16      0.18     0.00     0.51 1.56        7
## sd(CommAtt2M1)                 0.09      0.11     0.00     0.32 1.53        7
## sd(CommAtt3M2)                 0.12      0.15     0.00     0.43 1.53        7
## cor(Intercept,CommAtt2M1)      0.26      0.59    -0.83     0.99 1.50        7
## cor(Intercept,CommAtt3M2)     -0.38      0.55    -1.00     0.78 1.51        7
## cor(CommAtt2M1,CommAtt3M2)    -0.40      0.55    -1.00     0.77 1.51        7
##                            Tail_ESS
## sd(Intercept)                    39
## sd(CommAtt2M1)                   29
## sd(CommAtt3M2)                   29
## cor(Intercept,CommAtt2M1)        39
## cor(Intercept,CommAtt3M2)        39
## cor(CommAtt2M1,CommAtt3M2)       32
## 
## ~Participant (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.36      0.18     0.01     0.53 1.54        7
## sd(CommAtt2M1)                 0.23      0.12     0.00     0.34 1.53        7
## sd(CommAtt3M2)                 0.30      0.16     0.00     0.45 1.53        7
## cor(Intercept,CommAtt2M1)      0.71      0.49    -0.70     0.99 1.51        7
## cor(Intercept,CommAtt3M2)     -0.76      0.44    -1.00     0.57 1.52        7
## cor(CommAtt2M1,CommAtt3M2)    -0.78      0.43    -1.00     0.57 1.52        7
##                            Tail_ESS
## sd(Intercept)                    27
## sd(CommAtt2M1)                   31
## sd(CommAtt3M2)                   29
## cor(Intercept,CommAtt2M1)        29
## cor(Intercept,CommAtt3M2)        28
## cor(CommAtt2M1,CommAtt3M2)       28
## 
## ~TrialNumber_c (Number of levels: 54) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.02      0.01     0.00     0.05 1.00     1085     1283
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      2.76      0.07     2.61     2.90 1.00      361      650
## CommAtt2M1     1.62      0.08     1.46     1.78 1.00      414      698
## CommAtt3M2    -2.15      0.06    -2.26    -2.02 1.01      417      808
## Modality1     -0.87      0.20    -1.24    -0.49 1.02      129      123
## Modality2      0.45      0.11     0.23     0.66 1.02      190      235
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.56      0.01     0.54     0.57 1.00     3809     1448
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# not that there are troubles with convergence! (default priors)

plot(fit_eff6)

plot(conditional_effects(fit_eff6), points = TRUE)

# looks ok

pp_check(fit_eff6, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# nice!!! this is definitely improvement

pp_check(fit_eff6, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# even this looks quite nice

Model 7: adding our own priors

Now we will try to fit the same model but with our priors

# check what priors need to be specified
get_prior(Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c),
                prior = priors_eff7,
                data = final_data)
##                   prior     class       coef            group resp dpar nlpar
##                  (flat)         b                                            
##                  (flat)         b CommAtt2M1                                 
##                  (flat)         b CommAtt3M2                                 
##                  (flat)         b  Modality1                                 
##                  (flat)         b  Modality2                                 
##                  lkj(1)       cor                                            
##                  lkj(1)       cor                        Big5                
##                  lkj(1)       cor            Expressibility_z                
##                  lkj(1)       cor                 Familiarity                
##                  lkj(1)       cor                 Participant                
##  student_t(3, 2.6, 2.5) Intercept                                            
##    student_t(3, 0, 2.5)        sd                                            
##    student_t(3, 0, 2.5)        sd                        Big5                
##    student_t(3, 0, 2.5)        sd CommAtt2M1             Big5                
##    student_t(3, 0, 2.5)        sd CommAtt3M2             Big5                
##    student_t(3, 0, 2.5)        sd  Intercept             Big5                
##    student_t(3, 0, 2.5)        sd                     Concept                
##    student_t(3, 0, 2.5)        sd  Intercept          Concept                
##    student_t(3, 0, 2.5)        sd            Expressibility_z                
##    student_t(3, 0, 2.5)        sd CommAtt2M1 Expressibility_z                
##    student_t(3, 0, 2.5)        sd CommAtt3M2 Expressibility_z                
##    student_t(3, 0, 2.5)        sd  Intercept Expressibility_z                
##    student_t(3, 0, 2.5)        sd                 Familiarity                
##    student_t(3, 0, 2.5)        sd CommAtt2M1      Familiarity                
##    student_t(3, 0, 2.5)        sd CommAtt3M2      Familiarity                
##    student_t(3, 0, 2.5)        sd  Intercept      Familiarity                
##    student_t(3, 0, 2.5)        sd                 Participant                
##    student_t(3, 0, 2.5)        sd CommAtt2M1      Participant                
##    student_t(3, 0, 2.5)        sd CommAtt3M2      Participant                
##    student_t(3, 0, 2.5)        sd  Intercept      Participant                
##    student_t(3, 0, 2.5)        sd               TrialNumber_c                
##    student_t(3, 0, 2.5)        sd  Intercept    TrialNumber_c                
##    student_t(3, 0, 2.5)     sigma                                            
##  lb ub       source
##             default
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##             default
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##             default
##   0         default
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0    (vectorized)
##   0         default
# Define the priors
priors_eff7 <- c(
  prior(normal(2.5, 1), class = "Intercept"),                    # Prior for the intercept (baseline Effort)
  prior(normal(0, 1.5), class = "b", coef = "CommAtt2M1"),           # Prior for CommAtt level 2
  prior(normal(0, 1.5), class = "b", coef = "CommAtt3M2"),           # Prior for CommAtt level 3
  prior(normal(0, 0.5), class = "b", coef = "Modality1"),                # Prior for Modality level 1
  prior(normal(0, 0.5), class = "b", coef = "Modality2"),                # Prior for Modality level 2
  prior(normal(0, 0.5), class = "sd", group = "Participant"),      # Random effect for Participant (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "Concept"),          # Random effect for Concept (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "Familiarity"),     # Random effect for Familiarity (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "Big5"),            # Random effect for Big5 (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "Expressibility_z"), # Random effect for Expressibility (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "TrialNumber_c"),   # Random effect for Trial Number (Intercept)
  prior(exponential(1), class = "sigma"),                           # Prior for residuals
  prior(lkj(4), class = "cor")
)
# fit model
fit_eff7 <- brm(Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c), 
                data = final_data,
                prior = priors_eff7,
                cores = 4)

saveRDS(fit_eff7, here("06_Modeling", "models", "fit_eff7.rds"))
beep(5)
fit_eff7 <- readRDS(here("06_Modeling", "models", "fit_eff7.rds"))

# summary
summary(fit_eff7)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 1 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.15      0.18     0.00     0.49 1.55        7
## sd(CommAtt2M1)                 0.09      0.11     0.00     0.31 1.53        7
## sd(CommAtt3M2)                 0.11      0.15     0.00     0.40 1.53        7
## cor(Intercept,CommAtt2M1)      0.23      0.48    -0.59     0.95 1.53        7
## cor(Intercept,CommAtt3M2)     -0.27      0.48    -0.97     0.53 1.53        7
## cor(CommAtt2M1,CommAtt3M2)    -0.28      0.48    -0.98     0.52 1.53        7
##                            Tail_ESS
## sd(Intercept)                    29
## sd(CommAtt2M1)                   29
## sd(CommAtt3M2)                   28
## cor(Intercept,CommAtt2M1)        31
## cor(Intercept,CommAtt3M2)        36
## cor(CommAtt2M1,CommAtt3M2)       29
## 
## ~Concept (Number of levels: 21) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.05      0.03     0.00     0.12 1.01      549     1664
## 
## ~Expressibility_z (Number of levels: 63) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.41      0.04     0.33     0.50 1.01      488
## sd(CommAtt2M1)                 0.57      0.05     0.48     0.67 1.01     1153
## sd(CommAtt3M2)                 0.33      0.03     0.26     0.40 1.01     1652
## cor(Intercept,CommAtt2M1)      0.81      0.09     0.58     0.94 1.01      313
## cor(Intercept,CommAtt3M2)     -0.92      0.04    -0.98    -0.82 1.00     1036
## cor(CommAtt2M1,CommAtt3M2)    -0.73      0.07    -0.85    -0.57 1.00     1690
##                            Tail_ESS
## sd(Intercept)                  1110
## sd(CommAtt2M1)                 1819
## sd(CommAtt3M2)                 2925
## cor(Intercept,CommAtt2M1)       320
## cor(Intercept,CommAtt3M2)      2005
## cor(CommAtt2M1,CommAtt3M2)     2606
## 
## ~Familiarity (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.15      0.18     0.00     0.49 1.54        7
## sd(CommAtt2M1)                 0.09      0.11     0.00     0.31 1.53        7
## sd(CommAtt3M2)                 0.11      0.15     0.00     0.41 1.53        7
## cor(Intercept,CommAtt2M1)      0.23      0.47    -0.57     0.95 1.53        7
## cor(Intercept,CommAtt3M2)     -0.27      0.48    -0.97     0.54 1.53        7
## cor(CommAtt2M1,CommAtt3M2)    -0.28      0.48    -0.98     0.53 1.53        7
##                            Tail_ESS
## sd(Intercept)                    37
## sd(CommAtt2M1)                   30
## sd(CommAtt3M2)                   30
## cor(Intercept,CommAtt2M1)        31
## cor(Intercept,CommAtt3M2)        29
## cor(CommAtt2M1,CommAtt3M2)       30
## 
## ~Participant (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.25      0.21     0.00     0.50 1.74        6
## sd(CommAtt2M1)                 0.15      0.13     0.00     0.32 1.73        6
## sd(CommAtt3M2)                 0.20      0.17     0.00     0.42 1.73        6
## cor(Intercept,CommAtt2M1)      0.45      0.50    -0.50     0.95 1.73        6
## cor(Intercept,CommAtt3M2)     -0.49      0.50    -0.98     0.49 1.73        6
## cor(CommAtt2M1,CommAtt3M2)    -0.51      0.50    -0.98     0.47 1.73        6
##                            Tail_ESS
## sd(Intercept)                   168
## sd(CommAtt2M1)                  171
## sd(CommAtt3M2)                  170
## cor(Intercept,CommAtt2M1)       113
## cor(Intercept,CommAtt3M2)       133
## cor(CommAtt2M1,CommAtt3M2)      125
## 
## ~TrialNumber_c (Number of levels: 54) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.02      0.01     0.00     0.04 1.00     1390     1638
## 
## Regression Coefficients:
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      2.75      0.07     2.62     2.88 1.02      535      960
## CommAtt2M1     1.61      0.08     1.45     1.75 1.01      666      992
## CommAtt3M2    -2.13      0.06    -2.25    -2.03 1.01      664     1105
## Modality1     -0.64      0.19    -0.98    -0.24 1.01      298      281
## Modality2      0.33      0.11     0.10     0.54 1.01      385      510
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.56      0.01     0.54     0.57 1.00     7878     3132
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# not that there are troubles with convergence! (default priors)

plot(fit_eff7)

# looks scheisse, divergent chains
plot(conditional_effects(fit_eff7), points = TRUE)

# nicht so gut

pp_check(fit_eff7, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# stays the same

pp_check(fit_eff7, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# 

Note that adding varying slope is maybe not necessary for all the current variables, namely Big5, Expressibility, TrialNumber and Familiarity.

We need to have them in the model to avoid confounds (see DAG), but let’s take them one by one

  1. Familiarity

As it’s influencing CommAtt and Effort, it should be there as fixed effect to block the causal path.

We also expect different baselines of effort depending on how familiar two people are, so varying intercept is accurate for this

Do we expect that familiarity affects people differently? Perhaps not necessarily. However, for some people it might be that familiarity leads to bigger willingness to be more effortful, for some it means that they know each other so well that they don’t need to put so much effort into things (they are more conventional in a sense)

There should also potentially be an interaction between Familiarity and CommAtt - high familiarity might produce high effort in first attempt, but in the following attempts there will be less needed

so

  • Familiarity
  • (1 + CommAtt | Familiarity) or interaction
  1. Big5

As it’s influencing CommAtt and Effort, it should be there as fixed effect to block the causal path.

We also expect different baselines of effort depending on how, for example, open people are, so that’s why we go for varying intercept.

Do we expect that personality traits affect people differently? Probably not

so

  • Big5
  • (1 | Big5)
  1. Expressibility

As it’s influencing CommAtt and Effort, it should be there as fixed effect to block the causal path.

Do we expect different baselines of effort depending on how the concept is expressible? Yes, probably

Do we expect that expressibility could influence effort differently? Not sure

Model 8: full adapted

priors_eff8 <- c(
  prior(normal(2.5, 1), class = "Intercept"),                    # Prior for the intercept (baseline Effort)
  prior(normal(0, 1.5), class = "b", coef = "CommAtt2M1"),           # Prior for CommAtt level 2
  prior(normal(0, 1.5), class = "b", coef = "CommAtt3M2"),           # Prior for CommAtt level 3
  prior(normal(0, 0.5), class = "b", coef = "Modality1"),                # Prior for Modality level 1
  prior(normal(0, 0.5), class = "b", coef = "Modality2"),                # Prior for Modality level 2
  prior(normal(0, 0.5), class = "b", coef = "Familiarity"),
  prior(normal(0, 0.5), class = "b", coef = "Expressibility_z"),
  prior(normal(0, 0.5), class = "b", coef = "Big5"),
  prior(normal(0, 0.5), class = "sd", group = "Participant"),      # Random effect for Participant (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "Concept"),          # Random effect for Concept (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "Familiarity"),     # Random effect for Familiarity (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "Expressibility_z"), # Random effect for Expressibility (Intercept)
  prior(normal(0, 0.5), class = "sd", group = "TrialNumber_c"),   # Random effect for Trial Number (Intercept)
  prior(exponential(1), class = "sigma"),                           # Prior for residuals
  prior(lkj(4), class = "cor")
)

fit_eff8 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z +  (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 | Big5) + (1 | Expressibility_z) + (1 | TrialNumber_c), 
                data = final_data,
                prior = priors_eff8,
                cores = 4)

saveRDS(fit_eff8, here("06_Modeling", "models", "fit_eff8.rds"))
beep(5)
fit_eff8 <- readRDS(here("06_Modeling", "models", "fit_eff8.rds"))

# summary
summary(fit_eff8)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 | Big5) + (1 | Expressibility_z) + (1 | TrialNumber_c) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.01      0.01     0.00     0.04 1.00     1736     1639
## 
## ~Concept (Number of levels: 21) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.02      0.01     0.00     0.05 1.00     1140     1639
## 
## ~Expressibility_z (Number of levels: 63) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.05      0.02     0.01     0.08 1.01      582      601
## 
## ~Familiarity (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.02      0.01     0.00     0.05 1.06       53
## sd(CommAtt2M1)                 0.22      0.11     0.00     0.34 1.53        7
## sd(CommAtt3M2)                 0.29      0.15     0.01     0.44 1.53        7
## cor(Intercept,CommAtt2M1)      0.27      0.34    -0.45     0.81 1.13       23
## cor(Intercept,CommAtt3M2)     -0.25      0.34    -0.79     0.47 1.13       22
## cor(CommAtt2M1,CommAtt3M2)    -0.73      0.41    -0.98     0.37 1.53        7
##                            Tail_ESS
## sd(Intercept)                   775
## sd(CommAtt2M1)                   33
## sd(CommAtt3M2)                   32
## cor(Intercept,CommAtt2M1)       170
## cor(Intercept,CommAtt3M2)       168
## cor(CommAtt2M1,CommAtt3M2)       29
## 
## ~Participant (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.02      0.01     0.00     0.05 1.08       38
## sd(CommAtt2M1)                 0.09      0.11     0.00     0.32 1.53        7
## sd(CommAtt3M2)                 0.12      0.15     0.00     0.42 1.53        7
## cor(Intercept,CommAtt2M1)      0.12      0.35    -0.55     0.77 1.16       18
## cor(Intercept,CommAtt3M2)     -0.11      0.36    -0.77     0.56 1.17       18
## cor(CommAtt2M1,CommAtt3M2)    -0.29      0.47    -0.97     0.53 1.53        7
##                            Tail_ESS
## sd(Intercept)                    73
## sd(CommAtt2M1)                   34
## sd(CommAtt3M2)                   30
## cor(Intercept,CommAtt2M1)        54
## cor(Intercept,CommAtt3M2)        58
## cor(CommAtt2M1,CommAtt3M2)       29
## 
## ~TrialNumber_c (Number of levels: 54) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.03      0.01     0.00     0.06 1.00     1223     1532
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            2.72      0.01     2.69     2.75 1.00     3194     2949
## CommAtt2M1           1.60      0.03     1.54     1.67 1.00     1625     2534
## CommAtt3M2          -2.13      0.04    -2.21    -2.04 1.00     1843     2688
## Modality1           -0.55      0.03    -0.62    -0.48 1.00     2696     2881
## Modality2            0.30      0.03     0.24     0.36 1.00     2952     3221
## Big5                 1.15      0.04     1.07     1.22 1.01      833     1880
## Familiarity          1.17      0.04     1.09     1.25 1.01      758     1738
## Expressibility_z     0.39      0.01     0.37     0.41 1.00     3393     2397
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.62      0.01     0.61     0.64 1.00     5500     2915
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# not that there are troubles with convergence! (default priors)

plot(fit_eff8)

# looks scheisse, divergent chains
plot(conditional_effects(fit_eff8), points = TRUE)

# nicht so gut

pp_check(fit_eff8, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

# stays the same

pp_check(fit_eff8, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

# 

Let’s do the same model, but familiarity will not have varying slope

Model 9: full adapted 2

priors_eff9 <- c(
  # Prior for the intercept
  set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
  # Priors for the main effects using Student's t-distribution
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
  set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
  set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
  set_prior("gamma(2,1)", class = "sd", group = "TrialNumber_c"),

  set_prior("gamma(2, 1)", class = "sd", group = "Participant"), # let's try gamma because it limits the negative values a bit better
  set_prior("gamma(2, 1)", class = "sd", group = "Concept"),
  set_prior("gamma(2, 1)", class = "sd"),
  
  # Prior for residual standard deviation (sigma)
  #set_prior("normal(0, 2.5)", class = "sigma", lb = 0),
  set_prior("gamma(2, 1)", class = "sigma"),
  prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)


fit_eff9 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z +  (1 + CommAtt | Participant) + (1 | Concept)  + (1 | Expressibility_z) + (1 | TrialNumber_c), 
                data = final_data,
                prior = priors_eff9,
                cores = 4)

saveRDS(fit_eff9, here("06_Modeling", "models", "fit_eff9.rds"))
beep(5)

Ok now I am making a bit better sense of it

so we definitely have our main predictor CommAtt

We have other predictors that we believe moderate the effect - that is modality, familiarity, expressibility, big5

we want to have also varying intercept and slope for some

  • participants - because we believe that their baseline effort varies, and that the effect of commatt on eff might vary across them

  • same for concepts

we want to have varying intercept for TrialNumber because there might be variation between earlier and later performances (either learning, or fatigue)

Making sense of priors

We need to look a bit closer at our priors so that it’s reasonable

priors_eff9 <- c(
  set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
  set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
  set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
  
  set_prior("gamma(2,1)", class = "sd", group = "TrialNumber_c"),
  set_prior("gamma(2, 1)", class = "sd", group = "Participant"),
  set_prior("gamma(2, 1)", class = "sd", group = "Concept"),
  set_prior("gamma(2, 1)", class = "sd"),
  
  set_prior("gamma(2, 1)", class = "sigma"),
  
  set_prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)

fit_eff9_priors <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z +  (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c), 
                data = final_data,
                prior=priors_eff9,
                family = gaussian,#
                sample_prior = 'only',
                cores = 4)
## Compiling Stan program...
## Start sampling
pp_check(fit_eff9_priors, type = "dens_overlay", ndraws = 100)

pp_check(fit_eff9_priors,
         type = "stat",
         stat = "mean",
         prefix = "ppd") +
  coord_cartesian(xlim = c(-20, 20)) +
  ggtitle("Prior predictive distribution of means")
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(fit_eff9_priors,
         type = "stat",
         stat = "min",
         prefix = "ppd") +
  coord_cartesian(xlim = c(-50, 20)) +
  ggtitle("Prior predictive distribution of means")
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# this is still of, the minimum value should be 0


pp_check(fit_eff9_priors,
         type = "stat",
         stat = "max",
         prefix = "ppd") +
  coord_cartesian(xlim = c(-20, 50)) +
  ggtitle("Prior predictive distribution of means")
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Let’s try priors with exponential(1) for sd, whether it will fix the negative values

priors_eff_exp <- c(
  set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
  set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
  set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),

  set_prior("exponential(2)", class = "sd", group = "TrialNumber_c"),
  set_prior("exponential(2)", class = "sd", group = "Participant"),
  set_prior("exponential(2)", class = "sd", group = "Concept"),
  set_prior("exponential(2)", class = "sd"),
  set_prior("exponential(1)", class = "sigma"),
  
  set_prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)

fit_eff9_priors_exp <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z +  (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c), 
                data = final_data,
                prior=priors_eff_exp,
                family = gaussian,#
                sample_prior = 'only',
                cores = 4)
## Compiling Stan program...
## Start sampling
pp_check(fit_eff9_priors_exp, type = "dens_overlay", ndraws = 100)

pp_check(fit_eff9_priors_exp,
         type = "stat",
         stat = "mean",
         bins = 50,
         prefix = "ppd") +
  coord_cartesian(xlim = c(-20, 20)) +
  ggtitle("Prior predictive distribution of means")
## Using all posterior draws for ppc type 'stat' by default.

pp_check(fit_eff9_priors_exp,
         type = "stat",
         stat = "min",
         prefix = "ppd") +
  coord_cartesian(xlim = c(-50, 20)) +
  ggtitle("Prior predictive distribution of minimal values")
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# this is still of, the minimum value should be 0

pp_check(fit_eff9_priors_exp,
         type = "stat",
         stat = "max",
         prefix = "ppd") +
  coord_cartesian(xlim = c(-20, 50)) +
  ggtitle("Prior predictive distribution of maximal values")
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Let’s try to tighten the priors even more , especially the sds that are still pushing the posterior to negative values (which is not possible)

priors_eff_t <- c(
  set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
  set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
  set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
  set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
  set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
  set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
  
  set_prior("normal(0.5,0.1)", class = "sd", group = "TrialNumber_c"),
  set_prior("normal(0.5,0.1)", class = "sd", group = "Participant"),
  set_prior("normal(0.5,0.1)", class = "sd", group = "Concept"),
  set_prior("normal(1,0.1)", class = "sd"),
  
  set_prior("normal(0.5,0.25)", class = "sigma"),
  
  set_prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)

fit_eff9_priors_t <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z +  (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c), 
                data = final_data,
                prior=priors_eff_t,
                family = gaussian,#
                sample_prior = 'only',
                cores = 4)
## Compiling Stan program...
## Start sampling
pp_check(fit_eff9_priors_t, type = "dens_overlay", ndraws = 100)

pp_check(fit_eff9_priors_t,
         type = "stat",
         stat = "mean",
         bins = 50,
         prefix = "ppd") +
  coord_cartesian(xlim = c(-10, 10)) +
  ggtitle("Prior predictive distribution of means")
## Using all posterior draws for ppc type 'stat' by default.

# this look okay

pp_check(fit_eff9_priors_t,
         type = "stat",
         stat = "min",
         prefix = "ppd") +
  coord_cartesian(xlim = c(-15, 10)) +
  ggtitle("Prior predictive distribution of minimal values")
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# this is better but we still see some 0s

pp_check(fit_eff9_priors_t,
         type = "stat",
         stat = "max",
         prefix = "ppd") +
  coord_cartesian(xlim = c(-10, 15)) +
  ggtitle("Prior predictive distribution of maximal values")
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# this looks reasonable
Check out: https://paulbuerkner.com/brms/reference/set_prior.html

Now let’s exchange the priors for the tighter ones

Model 10: full with mildly informative priors

fit_eff10 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z +  (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c), 
                data = final_data,
                prior=priors_eff_t,
                family = gaussian,
                iter = 4000,
                cores = 4)

saveRDS(fit_eff10, here("06_Modeling", "models", "fit_eff10.rds"))
# warnings with ESS for correlation coefficients!
fit_eff10 <- readRDS(here("06_Modeling", "models", "fit_eff10.rds"))

# SUMMARY
summary(fit_eff10)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c) 
##    Data: final_data (Number of observations: 5055) 
##   Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.04      0.02     0.01     0.09 1.00     2037
## sd(CommAtt2M1)                 0.28      0.06     0.17     0.42 1.00     1824
## sd(CommAtt3M2)                 0.26      0.07     0.15     0.41 1.00     2139
## cor(Intercept,CommAtt2M1)      0.25      0.37    -0.53     0.85 1.01      625
## cor(Intercept,CommAtt3M2)     -0.20      0.38    -0.84     0.58 1.00      676
## cor(CommAtt2M1,CommAtt3M2)    -0.92      0.07    -0.99    -0.71 1.00     3026
##                            Tail_ESS
## sd(Intercept)                  1928
## sd(CommAtt2M1)                 3836
## sd(CommAtt3M2)                 3584
## cor(Intercept,CommAtt2M1)      1090
## cor(Intercept,CommAtt3M2)      1397
## cor(CommAtt2M1,CommAtt3M2)     4946
## 
## ~Participant (Number of levels: 120) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.06      0.02     0.03     0.09 1.01     1067
## sd(CommAtt2M1)                 0.32      0.03     0.27     0.38 1.00     3558
## sd(CommAtt3M2)                 0.41      0.04     0.34     0.49 1.00     2883
## cor(Intercept,CommAtt2M1)      0.82      0.13     0.47     0.97 1.01      467
## cor(Intercept,CommAtt3M2)     -0.82      0.13    -0.97    -0.46 1.02      468
## cor(CommAtt2M1,CommAtt3M2)    -0.97      0.02    -1.00    -0.93 1.00     4931
##                            Tail_ESS
## sd(Intercept)                   914
## sd(CommAtt2M1)                 5393
## sd(CommAtt3M2)                 4748
## cor(Intercept,CommAtt2M1)       357
## cor(Intercept,CommAtt3M2)       323
## cor(CommAtt2M1,CommAtt3M2)     6177
## 
## ~TrialNumber_c (Number of levels: 54) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.04      0.01     0.01     0.07 1.00     2168     2276
## 
## Regression Coefficients:
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            2.72      0.02     2.69     2.75 1.00     2992     3779
## CommAtt2M1           1.53      0.07     1.38     1.67 1.00     1816     3049
## CommAtt3M2          -2.05      0.08    -2.19    -1.89 1.00     2016     2732
## Modality1           -0.54      0.03    -0.60    -0.49 1.00     6604     6447
## Modality2            0.30      0.03     0.25     0.35 1.00     8200     6517
## Big5                 1.08      0.04     0.99     1.16 1.01     1662     2427
## Familiarity          1.09      0.05     1.00     1.18 1.01     1284     1397
## Expressibility_z     0.39      0.01     0.37     0.41 1.00     4057     5432
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.62      0.01     0.61     0.63 1.00    11144     5631
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# PLOT COEFFS
plot(fit_eff10)

plot(conditional_effects(fit_eff10), points = TRUE)

# PP CHECK
pp_check(fit_eff10, type = "dens_overlay", ndraws = 100)

# the fit looks quite ok

## PP CHECK WITH SUMMARY STATS
pp_check(fit_eff10,
         type = "stat",
         stat = "mean",
         bins = 50)
## Using all posterior draws for ppc type 'stat' by default.

pp_check(fit_eff10,
         type = "stat",
         stat = "min",
         bins = 50) # ok here we see it's looking for values where it should not
## Using all posterior draws for ppc type 'stat' by default.

pp_check(fit_eff10,
         type = "stat",
         stat = "max",
         bins = 50) # here it's also quite off
## Using all posterior draws for ppc type 'stat' by default.

## PP CHECK FOR RESIDUALS
pp_check(fit_eff10, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.

Where to next

What are potential interactions?

  • Modality*Expressibility (only highly expressive concepts might give modality the power to affect effort)
  • Big5*Modality (only in vocal modality traits could matter)
  • Big5 and/or Familiarity * CommAtt

Nesting

Dyad:Participant Modality:Concept Modality:Expressibility -

possibly change family to student , family = student,